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Vortex tilt response in a Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) vortex lattice is studied as a 
probe reflecting the spatial structure of this state. In quasi 2D materials under a parallel magnetic 
field, the tilt modulus E2 of the nodal planes in the FFLO state modulating along the field decreases 
as the paramagnetic effect is effectively enhanced, and this reduction of E2 in turn reduces the vortex 
tilt modulus. This reduction, more remarkable in higher fields or in more 2D-like systems, of vortex 
tilt modulus upon entering the FFLO state may be one origin of an anomalous anisotropic reduction 
of sound velocity detected in an ultrasound measurement in CeColn5. 
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I. INTRODUCTION 
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Recent accumulating experimental facts in heat capacity [1], thermal conductivity [2j], penetration depth 3], and 
NMR data 0| certainly indicate the presence of a new superconducting (SC) phase of CeCoIn5 at low T and under 
high magnetic fields parallel to the SC layers. This new phase, separated from the ordinary Abrikosov vortex lattice 
via a second order transition, is expected to be the FFLO vortex lattice with one-dimensional periodic modulation 
of the SC order parameter A parallel to the field H based on a consistency between these observations and a recent 
theory [B|, 0| on the characters of transitions between different phases. If the FFLO modulation is perpendicular 
to H, the mean field if ^-transition between the normal and the FFLO state is usually expected to be of second 
order just like the conventional result in the Pauli limit Q- However, this is incompatible with the discontinuous 
if C 2-transition [l], 0] in CeCoIn5. Spatial structures of a FFLO state may also be reflected in its elastic properties, 
and the tilt response of vortices should be sensitive to the direction of the periodic modulation. 

In this paper, we examine changes of vortex tilt modulus occurring through the transition between the FFLO vortex 
lattice with modulation parallel to H and the familiar Abrikosov lattice and show that, through a coupling between 
the vortices and the nodal planes accompanying the FFLO modulation , a measurable reduction of vortex tilt modulus 
may occur in such an FFLO state of uniaxially anisotropic superconductors in fields parallel to the SC layers. The 
present result may be relevant to the ultrasound experiment in CeCoIn5 [T3|, in which a monotonous reduction of 
sound velocity upon cooling through f?FFLO (T) was observed only for a sound mode accompanied by vortex tilts. 

First, a qualitatively expected feature of the tilt response in the FFLO state will be explained. The FFLO state of 
our interest is the so-called LO state in which A in equilibrium has a periodic modulation with a period 2ir/Q parallel 
to H || x and vanishes on periodic nodal planes lying in y-z plane (see Fig.l). The continuous FFLO to Abrikosov 
transition at -ffFFLo(T') is characterized by a vanishing of the Q 2 term in the mean field free energy. On the other 
hand, the vortex line tension 5C44 of the Abrikosov lattice, which is one part of its total tilt modulus 

B 2 

C44 = (5C 44 + — , (1) 

is defined from the gradient term in the fluctuation free energy for variations parallel to H, where B is the uniform flux 
density. Then, as in the conventional elastic softening in the ordinary solids, it is natural to expect H a tilt softening 
to occur on i?FFLo(T) with a cusp- like minimum of SC^ (see Fig. 2). Although such a behavior is generally expected 
for the so-called FF state with no nodal planes of the amplitude |A| and presumably also for other modulated vortex 
lattices with nodal planes parallel to the vortices, this picture is, when applied to the LO state with periodic nodal 
planes perpendicular to the vortices, justified only in the limited case where the nodal planes are never coupled with 
the vortices. In the LO state, a tilted nodal plane can carry the magnetic flux, and hence, a vortex tilt is induced by a 
small tilt of nodal planes since the number of vortices should remain unchanged for a small variation. This statement 
is represented in terms of tilt angles V and 9 n of the vortices and of the nodal planes (see Fig.l), respectively, by the 
clastic energy 

where A^(0) is the density of states per spin in the normal state, and the bracket ( ) denotes the spatial average. 
The dimensionless vortex line tension defined under fixed nodal planes is given by E±, E2 is the corresponding tilt 
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FIG. 1: Sketch representing tilted vortex lines (thick solid lines) and nodal planes (thin solid ones) in the FFLO state modulating 
along H of a quasi 2D system in a parallel field (H || x). Thick and thin dashed lines denote their positions in equilibrium. 
The SC layer and the anisotropy axis correspond to the x-y plane and the z-axis, respectively. 

modulus of the nodal planes, and a nonvanishing coefficient E3 of the coupling term is a consequence of the periodic 
modulation of the equilibrium A. Hereafter, the nonlocality arising from the long interaction range between the 
vortices is neglected in C44 [llj] . Then, the vortex line tension SC44 in the FFLO state stabilized by a positive E 2 
becomes 



<5C 44 = iV(0)T c 2 



E x -4 



En 



(3) 



Equation §5§ implies that 8C44 is reduced more drastically with decreasing E 2 . Such a decrease of E2 occurs due to 
an effective reduction of the orbital depairing because the modulation parallel to H is supported by the orbital effect 
of the magnetic field. In quasi 2D superconductors under a field parallel to the SC layers, the orbital depairing effect 
becomes less important in more 2D-like systems where the paramagnetic depairing is relatively more important. In 
the Pauli limit with no orbital depairing, the direction of modulation is spontaneously chosen, i.e., Ei is vanishingly 
small, as far as a Fermi surface (FS) anisotropy is negligible. Thus, in highly 2D- like systems and/or a case with a 
larger Maki parameter ojm, the FFLO state in H || x should show a softer tilt response, for displacements || y, as a 
consequence of a large fluctuation of nodal planes induced by an enhanced paramagnetic depairing. 

The physical argument given above implies that, in extremely 2D-like systems in the parallel fields, SC44 in the 
FFLO state may take a negative value and suggests a possibility that even a tilt instability of the FFLO state might 
occur. To sec to what degree the reduction of C44 is substantial in real systems, a consistent and microscopic derivation 
of E n (n = 1, 2, and 3) and the phase diagram will be performed in the remainder of this paper, and the ultrasound 
data [10( will be discussed based on the resulting tilt response. 

This paper is organized as follows. In sec. II, details of the model and our calculations performed to obtain the 
phase diagram and the tilt modulus are explained. In sec. Ill, typical examples of our numerical results following from 
the expressions derived in sec. II and Appendices are shown. In Appendix A, expressions following from a higher LL 
mode leading to a correction to H C 2 are given, and a theoretical background on the relation between the vortex tilt 
deformation and a LL mode of the order parameter is explained in details in Appendix B. 
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FIG. 2: Example of t (= T/T c ) dependence of L44, proportional to 8C44, in H || z which has a cusp-like minimum at the FFLO 
to Abrikosov transition temperature 0.288T C . 
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II. MODEL AND DESCRIPTION 



Our approach for describing the vortex lattices affected by the paramagnetic depairing takes the same route as the 
previous one [5| and starts with a BCS Hamiltonian H = Ho + Hj + Hint f° r quasi 2D systems, where 



Ho = dJ2 d 2 r ± (tf(r ± )y 



(-iV ± + eA) 2 
2m 



(4) 



H, 



and 



Wint — — 



(5) 



(6) 



with B^^kjJ = J2 P± ^p a j P-) a j (P+)> w here p± = px ± kj_/2. Here, j is the index numbering the SC layers, 
Px is the component of p parallel to the layers, A p is the normalized orbital part of the pairing-function which, in 
the case of d x 2_ y 2 -pairing, is written as y2(p%. — Py) in terms of the unit vector p parallel to the layers, and m is the 
effective mass of a quasi-particle. Further, ap,B = pB or —fiB is the Zeeman energy [12j . In discussing our calculation 
results, the strength of the paramagnetic effect will be measured by the Maki parameter au — \/21?orb(0)/ifp(0), 
where H mh (0) and H P (0) = irT c /(V2e~< E p,) ~ 1.2T c //i are the orbital and Pauli limiting fields at T = defined within 
the weak-coupling BCS model, respectively, where 7e = 0.577 is an Euler constant. Hereafter, the gauge field A will 
be assumed to consist only of the contribution of the uniform flux density B, i.e., we work in the type II limit with 
no spatial variation of flux density, because we are interested mainly in the field region near H C 2. 
We use the familiar quasi-classical approximation for the single-particle propagator 



Here, the Fourier transform of G £iCr (r) is given by 



/\ ie / tZs-A 



(7) 



G e ,<x(p) = ie + vfJ-B - e p 



(8) 



where e p = (p ± — p 2 F )/(2m) — Jcos(p 2 s), pf is the Fermi momentum in 2D case, and e„ denotes the Matsubara 
frequency 2nT(n + 1/2). Since we take account of the paramagnetic depairing suppressing the upper critical field in 
the mean field approximation H C 2(T), the use of the quasi-classical approximation, valid if pftb 3> 1, is safely valid, 
where tb — (2|e|£?) -1 / 2 is the magnetic length. 

Hereafter, in deriving an appropriate Ginzburg-Landau (GL) functional in H || x, the spatial variation of the SC 
order parameter field A in the out-of-plane (z-) direction is assumed to have a longer range than the interlayer spacing 
d. When the paramagnetic effect is absent, this continuum approximation is valid only if B <^ B CI = l/(2e~/d 2 ) [l3T |. 
Here, 7 (> 1) is the anisotropy defined from the ratio between the in-plane coherence length £0 — vp / [2i:T c ) and the 
out-of-plane coherence length £00 where vf is the Fermi velocity in 2D case. In the present case where the H C 2(T) is 
reduced via the paramagnetic depairing, this continuum approximation is safely valid if Hp(0) <C B CI , or 



7 < 



1 



2eH P d 2 



l-3£p (2D) 
d 2 a M 



(9) 



where a M = V2H^ ih (0)/f/p(0), and H^. h (0) is the orbital limiting field in 2D limit for fields perpendicular to 
the layers. The above inequality means that an increase of 7 competes with an enhanced paramagnetic depairing. 
In fact, eq.© implies that, as the paramegnetic depairing is stronger under a fixed anisotropy, the FFLO state just 
below H C 2{T) tends to enter not the nearly 2D region in B > B cr but the anisotropic 3D regime below B cl in which 
repeated structural transitions between the Josephson vortex lattices occur. Nevertheless, the transitions between 
the ordinary Josephson vortex lattices are not visible in most of quasi 2D materials, and the layer structure may be 
treated as an anisotropic 3D-like medium for most purposes as far as the intrinsic pinning effect of vortices does not 
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become essential 14[. Then, the difference e p+ — e p _ of the quasiparticle energy may be replaced by w • k even in 
layered systems. Consequently, the quadratic term of the GL free energy density is given by 



T 2 = - J d*rA*(r) \j- - K 2 {TL) 1 A(r), (10) 

where w is the velocity field on the FS, II = — iV + 2eA, and an appropriate gauge transformation has been performed 
to make the gauge held parallel to the z-direction. The operator K 2 is given by 



K 2 (IL) = ~£ / |A p | 2 G e , CT (p)G_ £ ,_ ff (-p + II) 

W Z^\l pi 2(ie n + afiH)- w-n/ FS 



N(0)J^ dpf(p)(\A p \ 2 exp(iT- 1 pwTl) 



FS 

(11) 



Here, ( )fs denotes the average over the Fermi surface (FS), t — T/T c , N(0) = N 2 (0)/d with the density of states 
A^ 2 (0) per spin in 2D case, 

f(p = . w 9 t r cos (J^LJL \ 12 
sinh(27rip) \ J c / 

and the identity 

i r°° 

-=j dpcM-pD) (13) 

was used in obtaining the last equality of eg. pip . Similarly, the 4-th order (quartic) term and the 6-th order one of 
the GL free energy density are written as 

Ti = h / d3 ^4n j -)A*(r 1 )A*(r 3 )A(r 2 )A(r 4 )| r .^ r , 

T& = h J rf3 ^6( n i) A *( r i) A *( r 3)A*(r 5 )A(r 2 ) A(r 4 ) A(r 6 )| r .^ r , (14) 



where 



K A = 



|E/ |A p | 4 G e , CT (p)G_ £i _ CT (-p + nDG_ £ ,_ (T (-p + n 2 )G £ , ff (p + n*-n 2 ) 



2^ ( o)TE( " iS8 ; ( ; ) ' Ap|4 \ , d5) 

f-'N dxd 2 d 3 / FS 



and 



K« = 



-?E / |A P | 6 G £ , ff (p)G_ £ ,_ ff (-p + n*)G_ £ ,_ ff (- P + n 6 ) 
x Ge,(T (p - nj + n 2 ) G- e ,-4- P + ni + n* - n 2 ) g £ , ct ( p - n 6 + n*) 

= T A,r (0)r y/ (- is g n ( £ ))l A p| 6 Ei.pcrm d t d l+1 d t+2 d z+3 \ 

^\ did 2 d 3 did 5 dQ (d\ + d 3 - d 2 )(d 3 + d 5 - di){d<s + d\ - d 6 ) / FS ' 

Here, X^perm did l+ \d i+2 d l+ z = did 2 d 3 d 4 + • • • + d 6 dxd 2 d 3 , dj — 2(k„ + cr/i.B) w IL, for an even j, and cZj = 
2(ie n + apB) — w • IT* for an odd j. Although the above expression of the 6-th order term is apparently different 
from the corresponding one in [5|, it can be numerically checked that both of them are the same as each other. By 
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using the identity (TT3|) again, we obtain 



K± = ^N(0) / n3 =1 ^ /(^^)(|A p | 4 exp 



,;=i 



— ( pi w • n 4 + p 2 w • n 2 + P3 w • rig 



FS 



#6 = ~^N(0) / nf =1 ^/£>)(|A p | 6 exp 



— Oiw • + p 2 w • n 6 



P3 w • (115 + n*- n 2 ) + P4 w • (n* + - n 4 ) + p 5 w • (n; + nj - n 6 ) ) 



(17) 



FS 



Hereafter, the order parameter field A will be decomposed into the LLs. When the FFLO state modulating along 
H is formed, no additional spatial variation is induced in the y-z plane perpendicular to H || x except that due to the 
vortex structure, and hence, the y and z dependences of A of the FFLO state in equilibrium is, as in the Abrikosov 
state, well described in the lowest (n = 0) LL. Nevertheless, in H || x parallel to the layers, the anisotropy 7 in the y-z 
plane between the coherence lengths needs to be determined to consistently define creation and annihilation operators 
II± for the LLs. To determine 7, we follow its derivation in the conventional GL region and focus on the quadratic 
term in II and II* which is proportional to 



((WyTLy 



w z TL z ) 



)fs = 7 [r^ + h^U v i7-^IL)( 7 ^n 



-1/2 T 



1/2 ^ + i7- 1/2 n z )], 



(18) 



where 7 = 



))/( w z)- Here and below, we have chosen the gauge A = —Byz leading to n_II + — n + II_ = 1 where 



n± = t\b(7 1 / 2 IT Z ± i7 1 / 2 n y )/i/2 are the creation and annihilation operators of the LLs satisfying H^fo(y, z) = 0, 
and (fn{y, z) is a basis function in the n-th LL. According to TLq defined above, the velocity field on the FS is given 

by 



t>, k z ) = Vp(l — J(l — cos(pzd))) 1 / 2 (cos((/))x + sm(4>)y) + Jdsm(p z d)z 



(19) 



where J = 2mJ /p 2 F . In this case, we have 7 = 2\ 1 — J /(nJ). For this w, the averaging over the FS is defined by 



(M) 



FS 



# r d( Pz d) 



2lT 



2tt 



M. 



(20) 



Next, the operation (p n {y,z) = exp(ipT c x w • Yl)cp n (y, z), necessary to make the expressions of Ti n tractable for 
numerical calculations, will be examined. This is most easily accomplished in terms of the corresponding coherent 
state exp(-H72)]T n 

>o °" ni Pn/ Vn\. Examining the action of exp(iT c 1 w • n) to this coherent state, we obtain 



(p n (y,z) = exp( i— w II )<p n (y,z) 



A/o_ 



—nr, exp — (p 



i/*r) p(p-p) + 



d{pp) 



exp — \pz 



(y + V2pp+pf 



for n < 2, where y = 7 1 ^ 2 y, z — 7 1 ^ 2 z, 



7 1/2 w y + i7 1/2 w z 
V2r B T c ' 



and the corresponding n — LL state is 



fo(y,z) = Af cxp( -ipz 



{y+pf 



(21) 



(22) 



(23) 



At this stage, it is straightforward to study properties in equilibrium of an FFLO vortex state. First, as already 
mentioned, we assume the FFLO state in equilibrium to be described in the lowest (n = 0) LL where no nodal points 
or lines except the field-induced vortices of |A| are present in the y-z plane perpendicular to H. Instead, nodal planes 
perpendicular to H are periodically formed. If the Q-dependence of the free energy is incorporated only from the 
quadratic term J~i , the transition line between the LO state and the Abrikosov state is the same as that between the 
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FF and Abrikosov states, where Q is the wavenumber corresponding to the period between the FFLO nodal planes. 
However, once the Q dependence of the free energy from the higher order terms, Tk {k > 4), is considered, as pointed 
out elsewhere [l6[ in the H || z case, the transition temperature between the LO and Abrikosov states is higher than 
that between the FF and Abrikosov states. For this reason, we will not consider a possibility of appearance of the FF 
state. Hereafter, the equilibrium order parameter A e in the FFLO state is assumed to take the form 



A e (r) = V2T c a e * A (y, z) cos{Qx) 



(24) 



with the normalization condition (|A e 



/sp 



= T 2 where 



/sp 



implies the spatial average, and \1/a is the Abrikosov 



solution of the vortex lattice in n = LL. Then, in equilibrium, the mean field free energy density of the FFLO state 
takes the form 



N(0)T? 



a (q)a 2 e 
c (0) («e) 



6 6 



.(4) 



(25) 



Hereafter, the dimensionless wavenumber of the FFLO modulation will be defined as 

M 2 



q 



'TBI 



(26) 



The coefficients in eq. ([25|) are given by 



Vi{q) 



a o (0) + a o V 



(4) 4 

% Q , 



= V 4 (0)-V^q 2 + V 4 



(4)„4 



(27) 



Hereafter, the (/ dependence of Vq will be neglected. This simplification is in part based on the fact that it has 
been verified [1 61 ] that the q dependences of Vq are unimportant even quantitatively for the position of the FFLO 
to Abrikosov transition in H || z. Further, to study systematically possible phase diagrams including FFLO states, 
inclusion of g-dependences in higher order terms requires a difficult numerical task. On the other hand, if even the q 
dependence of V4 is neglected, the B-T region in which the FFLO state can appear is highly overestimated, and, as is 
seen in sec. Ill, a fictitious tilt instability of the FFLO state occurs. Therefore, for the practical purposes, neglecting 
q dependences of Vq and keeping the corresponding ones of Vi(q) is a convincing approach. Of course, when using 
eq. (P5)l . it is necessary to verify the conditions Vq > and > which justify the use of T e truncated at the 
0(|A| 6 ) and 0(g 4 ) terms. 

The onset temperature To at which the mean field i/ C 2-transition becomes discontinuous is given as the position at 
which Vi(q) becomes negative upon cooling while Vq > 0, and a second order transition line Hfflo{T) is determined 
as the line on which c^ 2 \a c ) in B < H C 2 becomes negative on cooling, while c^(a ) > 0, . Then, by minimizing T e 
with respect both to q and a e , a 2 is determined by 



a 2 e (q) = 

while q = if a (2) - V 4 (2) (a e (q)) 2 /2 > 0, and 



-VM + y / (^4(g)) 2 -4a (g)1/ 6 
2Vb 



-a^+v} 2) (a e (q)) 2 /2 



2(-a^ + (a e (q)) 2 Vr/2)- 

otherwise. Below To, the discontinuous iJ C 2-transition (i.e., first order mean field SC transition) occurs when 

, s 3 (^ 4 (g)) 2 
ao((?)= 16^^ 



(4) 



(28) 



(29) 



(30) 



irrespective of the minimized value of q. 

By applying eq. ([2"4")) to Tj, it is straightforward to derive the coefficients in eq. p5"|) if eq. p^)) is repeatedly used. 
Consequently, they are given by 



oo(0) = hn(h)+ f d P 
1 Jo 



1 ( ^ey 

— exp ' 

P 



^,-/(p)(|A p | 2 exp 



\P?P 2 



FSJ 



7 



,(2) 



dpf{p)p\ (Re/i) cot^|A p |"exp - 



IpIV 



FS 



1 

6 Jo 



dpf(p)p\ (Rep) 4 cot 4 0|A p | 2 exp 



Vi(0) = 



(2) 



poo / 3 

3/ If? :<^/ E 

/• oc / 3 

3/ ir E 



Pj 



Pj 



FS 
i?14 



COs(/4) 



FS 



E^-E^: 



i+j 



P1P3 



(Rep) 2 cot 2 (/>|A p | 4 exp 



3=1 

1 

r 3 



i?24 + -Rl4 ) ) C0S(I 4 ) 



FS 



1 f°° t \ r 4 

^0 v j=1 n j=1 ^ 3 



(Rep) 4 cot 4 0|A p | 4 ex P -- \--R 



3=1 
1/1 



COs(/4) 



FS 



-15 y n j 5 =1 dp J /^ Pfe ^|A p | 6 exp^i(i? 16 + i? 2 
where ft = B/H^P\o), and 



cos(/ 6 ) 



(31) 



FS 



i?14 
h 

R26 



ex 
ei 
e3 
e 4 

es 
e6 

eg 
eg 



H 2 (E^ + P2(p 3 +Pl)), 
3 = 1 

Re(p 2 )(P2 + (P3-pi) 2 ), 



Mm 2 )/ j 



\p\ \ ex + e 2 + e 3 + -e 4 e 5 



3 

el -el 



Rep ^ex + e 2 + e 3 
Im(p 2 ) 

^ (ex + e 2 - e 3 H — 

(p3 + ,05) 2 + (A3+P4) 2 , 

(Pi + P4+ P5f, 

P3 + P4 + (Pa - P5) 2 , 
px + 2(p 3 + p 4 + /05), 

Pi - P3 - Pi- Ph, 

(P4 - Ps) 2 + (Pi + P5 ~ P3) 2 , 

(Pi +P4 - P3) 2 , 

(/03-/0 4 ) 2 + (P2+P3-P5) 2 , 

(Pa + P4 - P5) 2 - 



-A ~T~ <^5 — , 

g(e6 + e 7 + e 8 + e 9 ) 



e 9 - e 6 - e 7 )) 



(32) 



In obtaining ao(0), the interaction strength \g\ has been eliminated under such a condition that, in H perpendicular 
to the layers, the operation — A" 2 at T = and in the absence of the paramagnetic effect vanishes at H^F'(0). 

It should be noted that, for any FS with anisotropy in the y-z plane, the n = 2 LL mode couples to the n = LL mode 
of A in high fields, and hence, that the expressions of GL coefficients given above are, strictly speaking, insufficient. 
This coupling inevitably occurs except in the conventional GL region valid in lower fields, where the gradient terms 
are kept only up to 0(II 2 ), and the y and z dependences of expressions were isotropized in determining 7 (see the 
description around eq. (fT5)) ). Expressions of coefficients in the GL quadratic term related to the n = 2 LL modes are 
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given in Appendix A. For the anisotropy values (7 < 5.5) used in our numerical calculations, however, this coupling 
was quantitatively negligible, and thus, the coefficients in ea. (|31[) are used in obtaining numerical results in sec. III. 

Now, let us turn to explaining how to describe tilt deformations in the FFLO state with nodal planes perpendicular 
to H. In the present parallel field configuration, the in-plane vortex tilt unaffected by the intrinsic pinning effect of 
the layering is expressed as a x-dependent vortex displacement u = uy related to a vortex flow in the y-direction. 
In a vortex lattice in equilibrium described in the n = LL, such a tilt deformation accompanied by a vortex flow 
is, consistently with the vanishing of the corresponding static superfluid rigidity, expressed as a fluctuation in the 
next lowest (n — 1) LL of A. [17] Such a relation between the n = 1 LL mode and the vortex displacement will 
be reviewed in Appendix B within the conventional GL region valid in low fields and near T c . Through the analysis 
shown there and in Ref . [13] , it is convincingly understood that the following points are valid beyond the conventional 
GL region. First, in examining the elastic deformations of the Abrikosov lattice, the energy gap between the n = 1 
and n = LLs is lost due to the magnetic screening, i.e., a gauge field fluctuation coupling to the vortex motion, 
and this disappearance of the mass gap is equivalent to the vanishing of the static superfluid rigidity, T s ^ = 0, for a 
phase twist perpendicular to H. The resulting main term of the vortex tilt modulus is the magnetic energy B 2 /(An) 
which is insensitive to the details of the SC state. Clearly, this result that the main term of C44 becomes insensitive 
to the details of the SC state as a consequence of the vanishing of T s ± holds true in the FFLO state modulating 
along H. On the other hand, according to the results in the conventional GL case in Appendix B, the remaining 
term SC44 = C44 — B 2 /(An) arises directly from the gradient (d x ) term of the resulting GL action regardless of the 
magnetic screening. That is, as far as focusing on 5C44, we can work in type II limit with no fluctuation of the gauge 
field. Then, consistently with eq. (|24|) . the SC order parameter field with tilt deformations of the vortices and of the 
nodal planes included should take the form 

A = A e (y, z) + SA(r) = V2a e T c ((p (y, z) + Sai(x)(fi(y, z)) cos(Qx + T c s{x, y)/v F ), (33) 

where s(x, y) is the displacement of the nodal planes, and the amplitude Sai(x) of the n = 1 LL fluctuation is identified 
with the vortex displacement u = uy parallel to the layers in the manner 

6ai(x) = ^-, (34) 



as explained in Appendix B. If s(x,y) = 0, SC44 is obtained as the coefficient of (d x u) 2 term. 

Hereafter, the elastic constants Ej (j — 1, 2, and 3) introduced in sec.l will be expressed in the manner 



E 1 = 2n 2 ^al(L 1 + ^-L 1 ). 

,2 



E 2 = a 2 e [L 2 + ^-L 2 y 

E * = ^^^ a e( L * + Y l3 )- (35) 

To first obtain the contributions Lj (j = 1, 2, and 3) from the quadratic GL term !F 2 , let us consider the following 
quantity 

(0|n) = 2(Lp* Q {y, z) cos{Qx + T c s{x, y)/v F ) expA^wII) 

x Sa n (x)(p n (y, z)cos(Qx + T c s(x,y)/v F )\ \ (36) 

/ FS/ sp 

(n = or 1) appearing commonly in Tj (j = 2, 4, and 6), where Sao = 1 is assumed. First, using the identity 
exp(^4 + B) = exp(A) exp(B) exp(— (AB — BA)/2), valid when AB — BA is a constant, the expression 

exp( -^pw • II )Sa n (x)ip n (y, z) cos(Qx + T c s(x,y)/v F ) (37) 



will be written as 



Sa n (x + T c 1 w x p)cos(Q(x + T C 1 pw x ) + T c s(x + T c 1 pw x ,y + T c 1 pw y )) exp ^-p( w v u y + w z ll z )^(p n (y, z). (38) 
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Then, the average on x over the scale 2ir/Q will be performed prior to all of the spatial averages by assuming a slow 
variation of s(x) in x, and thus, 2cos 2 (Qa; + s(x,y)) may be replaced by unity. Further, using eqs. ([37|) and (|38f and 
keeping only terms remaining finite after the momentum average on the FS, we find 



and 



(0|0) = (cos(pQw x /T c ) (cos(pwyd y s/v F ) cos{pw x d x s/v F ) exp(-p 2 |p| 2 /2)) sp ) FS , 

pM 2 



(39) 



p-^-d x 8a 1 {x) 

-L c 



w 

VF 



y pd y s(x,y) ) sin( pQ^- ) (-p/i*) exp 



sp/ FS 



2^ 2 p 3 

7 



CO Wy 



r B v F J v F 



->x 



w 



sin( p—Q ) exp 



p 2 |p| 2 



(d x ud y s} S p. 



(40) 



l-'S 



The latter is valid up to the harmonic order in s and u. Then, Z^'s are easily obtained in terms of the above expressions 
and are given by 



L 3 



dp [\A p \ 2 f( Po ) 
d P J |Ap| 2 /( Po ) 



w x p 
v F 

WyPO 
VF 



(l- P ^|/i| 2 )ex P - 



P§lM| 2 



cos(0) 



(41) 



FS 



exp 



dpa{\Ap\ f{pa)Re{p) p — exp 



cos(0) 



FS 



sin(0) 



FS 



where 



cos(n) = cos^\/2 qp n Re(p,)cot</> 
sin(n) = sinl V2qp n Re(/i)cot</> 



(42) 



It might be natural to discuss the elastic deformation of the vortex lattice based on these Lj's without including the 
contributions from JF4. However, L 2 itself is found not to lead to a qualitatively reasonable result of E 2 in the FFLO 
state: As is shown later in Fig. 4, the resulting tilt rigidity L 2 of the nodal planes often becomes negative. This result, 
suggestive of an instability of the FFLO state modulating along H, is an artifact due to the neglect of contributions to 
E 2 from the higher order terms, J-4, of the GL free energy. Hereafter, consistently with the neglect of g-dependences 
in Vq in eq. ([23|) . E n will be expressed, as in eq. ([55|) . as the sum of the contributions of T 2 and J-4 terms of the GL free 
energy. Derivation of the contributions Lj to the elastic moduli from T± is lengthy but straightforward once using 
the expressions (|37| and (j38|) , and they are expressed by 



U 3 3=1 dp J (px+p 2 ) 2 f(Y / P 3 )(\A p \ i (^)Y3n3 =1 cos(j) + i^(-l)^sin(z)sin(j)cos(6 



-1 + - [ (3px +p 2 - p 3 )(3p 2 + Px+ p 3 )(|p| 2 + Re( M 2 )) + [(3/9i - Pa - P 3 )(3p 2 -pa- Pi) + 4(pi - p 2 ) 2 



+ 4p 3 (p! + p 2 )](|p| 2 - MPV) ] ) cos(h) + -Im(p) 2 ((p 2 + p 3 ) 2 - pf) sin(/ 4 ) 



exp 



1 f 1 



-R 



24 



R 



14 



FS 



l 2 = f n^dft./fcp/j/iA/f^)^ 



+ (-l) 4+:, cos(6 - % -j) ( p l p ] cos(i)cos(j) - -{pe-i-j) sin(i)sin(j) 



exp 



1 f 1 



i?24 + R 



14 



i-j) 
cos(/ 4 ) 



Pj pi 



Jo \ j=1 
— 2(— l) 4+3 cos(6 — i — j) Pi cos(i)sin(j) 



4 w x w y 



Pl vl 



E 



2 v 2 

Pe-i-3 sin(6 - i - j)( 3cos(z)cos(j) + (-l) i+ - 7 sin(?)sin(j) ) 



FS 



Re(p) (3pi - p 3 + p 2 ) cos(I 4 ) + Im(p) ftj si 11 ^)^ 
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x exp^-i^-ii? 24 + i?i4^ • (43) 

Actual numerical calculations of tilt moduli are performed according to eq. (|35p in terms of Lj and Lj (J = 1, 2, and 
3) given above. 

Before ending this section, we point out a couple of essential features appearing in the coefficients of the GL free 
energy and the elastic constants derived above. First of all, noting that the contributions of the orbital depairing 
appear as \^\p n or Re(^)p„ everywhere, the effective strength of the paramagnetic depairing is cum (h 1 ) 1 / 2 . Hence, an 
increase of the flux density B or of the anisotropy 7 enhances the paramagnetic depairing effects. Further, consistently 
with this discussion, the period of the FFLO modulation is scaled by not the coherence length but the magnetic length 
tbI 1 / 2 (see cq. (j2"r?)) ), implying that the period of the modulation decreases with increasing B. Although this is not 
surprising because the paramagnetic effect is enhanced with increasing B, one should note that rs does not arise in 
any approach in the Pauli limit with no vortices included. 



III. NUMERICAL RESULTS AND DISCUSSIONS 



In this section, typical numerical results following from the expressions of the coefficients, cq. ([3"Tj) . and the elastic 
constants, eq. (|35j) . are presented, and their relevance to an available experimental result will be discussed. 

In the ensuing numerical results, the values olm = 10.65 and 7 = 4.5 were used for H || x, otherwise stated. Further, 
for comparison, clastic constants in perpendicular fields, H || z, were also examined in terms of the values olm — 4.95 
and 7 = 2.85. It is straightforward to, with an appropriate replacement of FS and w, obtain the corresponding 
expressions to eqs.(f3~Tj) and (f35|) in perpendicular fields which were used elsewhere [H, [lfjj]- The orbital- limiting field 
-fforb(O) in each field configuration was estimated numerically from ao(0) — with olm = in T — ► 0. 

In determining the phase diagram Fig. 3(a), the onset temperature To of the discontinuous iJ C 2-transition was 
determined as the position at which V±{q) changes the sign while verifying Ve(q) > 0, and i?FFLO (T) is determined, 
when cW(a e ) > 0, as the line on which c^(a c ) in H < H C 2 changes the sign. The above-mentioned conditions on 
the sign of Vq and were satisfied in all of the resulting numerical data in H || x. As reported elsewhere [lij ]. 
the second order transition on Hpflo(T) occurs at lower temperatures than To and decreases upon cooling, since, 
as mentioned at the end of sec. II, the paramagnetic depairing is stronger in higher fields and at lower temperatures. 
We note that, in the case of Fig. 4 (a), the FFLO state modulating along H is overcome, in t < 0.25, by another 
FFLO-like vortex state with a modulation perpendicular to H and formed in the next (n = 1) LL However, 
this another FFLO-like state has no periodic modulation parallel to H, implying that no specific feature is expected 
in tilt deformations. For this reason, we focus here on the higher temperature range in which the n = 1 LL state does 
not occur. 

The main result in this paper is seen in Fig. 3 (b): As the curve of L^{T) in h = 0.5 shows, where 




FIG. 3: (a) Example of the h-t mean field phase diagram in H || x obtained numerically, where each thick (thin) solid curve is 
the discontinuous (2nd order) mean field transition curve. The low temperature region in t < 0.25 where another FFLO-like 
vortex lattice [1, [l8j] described by the n — 1 LL modes of A e occurs is not shown here, (b) The corresponding numerical data 
of L44, defined in the text, in clean limit. The lower (upper) solid curve denotes Li4,i(t) in H || x for h — 0.5 (h = 0.485), while 
the dashed curve is that in H || M given in Fig. 2 and follows from the dotted curve in Fig.4. 
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-0. 004 



0.3 t 34 



0.26 



FIG. 4: Numerical results of L2 in H || x for h — 0.485 (upper solid curve) and h — 0.5 (lower solid one) and in H || z for 
h — 0.24 (dotted one) obtained by including contributions from the GL-quartic term to the description of the FFLO state. The 
dashed line showing L2 < follows from ea. (|42p with 7 = 3.65. 

6C44 proportional to L44 (see eq.®) can become negative in the FFLO state and thus, leads to a reduction of C44. As 
argued in sec. I, this 5C44-reduction is more remarkable in higher B and at lower T, i.e., as the paramagnetic depairing 
is stronger. 

In Fig. 4, results of the nodal plane's tilt modulus E2 in various situations are shown. When L2 is neglected in E2, 
and only the contribution L2 from T2 is kept, we often see, as in the dashed curve, negative E2 values. It implies that 
theoretical approaches based on an evaluation of an inhomogenuity of the order parameter field A(r) only from the 
quadratic GL term T2 [HI cannot describe the stability of the FFLO state properly as far as the normal to FFLO 
transition is discontinuous in the mean field approximation. In contrast, once L 2 is included in E2, the results of 
E 2 we obtain remain positive, as is physically required for a stable FFLO state modulating along H, although they 
significantly decrease upon cooling especially in higher fields, reflecting the "softening" of nodal planes induced by 
the paramagnetic depairing (see sec. I). Note that the decrease of E2 is not unlimited according to our microscopic 
calculation, suggesting that a genuine instability of the FFLO state modulating along H does not necessarily occur. 
In any case, the primary origin of the negative SC44 in the FFLO state is clearly this reduction of E2 upon cooling 
which becomes more remarkable in situations affected by a stronger paramagnetic depairing realizable in H || x. In 
this way, the physical consideration given in sec. I is supported by the microscopic derivation of the elastic moduli in 
sec. II. 

The reduction of C44 is estimated based on 



The expression in the bracket [ ] is evaluated as 3.5 in terms of the Li-result known within the weak-coupling BCS 
model. Based on Fig. 3 (b), R can reach —0.5 in the temperature range just below Hfflo(T). Imagining CeCoIn5 
in several tesla and, as a rough estimation, taking the values A(0) ~ 10 3 (A) and tb — 10 2 (A), the resulting —SC44 
values in the FFLO state near TJfflo^) are two order of magnitude smaller than B 2 /(4ir). This small value of 
reduction of 8C 44 /C 44 will partly become larger by taking account of the dispersive main term [17[ of C44 (the k 
dependences of the B 2 /4n term). 

The present result may be relevant to the observation in the ultrasound measurements [Io| . where the normalized 
sound velocity v s has shown a reduction upon cooling through -Hfflo for the displacements of the underlying crystal 
perpendicular to H (Lorentz mode), while no signature of a comparable magnitude has been seen for displacements 
parallel to H (nonLorentz mode). This remarkable anisotropy of phenomena is an evidence of a structural change 
on Hfflo(T) of the vortex state. Since, strictly speaking, a pinning of the nodal planes due to the crystal lattice is 
present in the latter, this fact implies that the pinning of the nodal planes is quantitatively negligible. In contrast, 
in the Lorentz mode, not the nodal planes but the vortices couple to the crystal displacement, and consequently, 
the observed reduction of v s in this mode in entering the FFLO state implies some reduction of C44 and/or of the 
vortex pinning strength. First, the overall temperature variation of v s surviving at low enough t likely reflects that 
of the order parameter amplitude a carried by the pinning strength (see Fig. 3 in Ref . (Io| ) . In addition, the data of 




2 • 



(45) 



where A(0) is the London penetration depth at T = 0, and 




(46) 
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quantities measuring a e (T) 0,0, 2(3, 21 1 show a reduction of a c (T), due to an increase of q (lij| . upon entering the 
FFLO state compared to the extrapolation of a e in the Abrikosov state to lower temperatures. Thus, the relative 
reduction of a e and hence of the pinning strength may be one origin of the ^-reduction in entering the FFLO state. 
On the other hand, the reduction of C44, obtained in the present work, upon entering the FFLO state, also leads to a 
reduction of v s in the Lorentz mode [10(. The feature seen in Fig. 3 that the C44-reduction is more remarkable in higher 
fields is consistent with the field dependence of the observed anomalies [l(J • Although the reduction of C44 estimated 
above seems to be too small to quantitatively explain the observed reduction of v s , we expect this discrepancy to be 
partly resolved by going beyond the GL expansion in A because a 2 has been presumably underestimated in the GL 
analysis of the present case with a discontinuous ff C 2-transition. 

For comparison, SC44 in H || z was also examined by adding a small noncylindrical portion with large \w z /w x \, 
stabilizing 0, [l8| an FFLO in H |j z, to the cylindrical FS. Although, as given in Fig. 2, it shows the familiar type of 
softening behavior of <5C44 at -HfflOj this cusp- like feature is not visible on the scale of Fig. 4 (see the dotted curve), 
and, as shown there, the temperature variations of SC44 and L 2 are quite weak around Hppho ■ Thus, changes of the 
tilt response through -ffFFLO in H || z are negligible. The presence in H || x and the absence in H || z of a clear 
peak effect in magnetization data of CeCoIn5 near -Hfflo [20| might be related to the corresponding difference in 
C44 mentioned above, because a peak effect occurs reflecting a notable change of a vortex elastic modulus [22j . 

Throughout this paper, we have focused on the case in which vortex lattices in equilibrium are described in the 
lowest LL. As shown elsewhere, [1, [H( a higher LL vortex lattice tends to occur in "clean limit", which is defined as 
the case with infinitely long quasiparticle's (QP's) mean free path, and in the case with large ajvf- In fact, in Fig. 3 
(a), the n = 1 LL vortex lattice with additional nodal lines in the plane perpendicular to H occurs in much lower 
temperatures than the range shown there. Such an appearance of higher LL states due to a large au under a strictly 
parallel field to the layers is closely related to a different issue [24| of transitions between 2D vortex lattices in the 
large 7 limit induced by a tilt of the applied field from the parallel field configuration. However, the n = 1 LL state 
in the present case is easily pushed down to T = due to a finite but nevertheless quite long QP's mean free path 
and is expected not to occur at measurable temperatures in CeCoIn5 In relation to this, it will be valuable to 
point out that the 2D higher LL vortex lattices due to the tilt of the applied field [24| may be expected only when 
the opposite relation to eq.(0), i.e., 



„ l-3£o (2D) 

7> -tf~ a M > 



(47) 



is satisfied. In fact, the neglect [24( of the orbital depairing effect in the case with a strictly parallel field is equivalent 
to assuming the absence of structural transitions between different Josephson vortex lattices, and this assumption is 
justified only in higher fields than l/(2ejd 2 ) [13j. For CeCoIn5 with a weak anisotropy [6(, eq. (|47f is never satisfied. 

In conclusion, the vortex tilt modulus C44 in the FFLO state modulating along H may be reduced due to tilts of the 
nodal planes. This C44-reduction should be remarkable especially in quasi 2D materials with a strong paramagnetic 
depairing in the parallel field, in which the nodal planes are fixed only weakly by the field direction, and may be an 
origin of the reduction of sound velocity in the Lorentz mode upon entering the FFLO phase observed in CeCoIn5. 

The author is grateful to Y. Matsuda and C. J. van der Beek for discussions on this issue. 

APPENDIX A 

Strictly speaking, ao(0) and other coefficients in eq. (|31[) . described within the n = LL, are affected by a mixing 
of the n — 2 LL modes in expressing the equilibrium order parameter A e . For instance, ao(0) should be replaced by 
a (0) - (a 2 o) 2 /a 2 (0), where 



a 2 (0) = Jln(/i)+ f dp 
* Jo 



O20 



1 



-(irZap/r B ) 



f(p)(\A p \' 1-2/l/z 



e 

_> 2 _> , / imIV 



exp 



ImIV 



dpf(p)( \A p \ 2 p 2 p 2 exp 



FS 



FS 



(48) 
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APPENDIX B 



We work here in the conventional GL free energy 
~T-T, 



•Ft GL — 



d 3 r 



J2 + 2eA,)#| 2 + J- 2 \Hd z + 2eA z )V\ 2 



^|4 + i. (curMA) 2 



(49) 



where b > 0, and A = Byz + SA(x). Here, the order parameter <3> is assumed to be the sum of the Abrikosov solution 
*A = a e * (0) in the n = LL and its excitation a e (5ai(x)* (1) in the n = 1 LL, where = II" /\fn\. Noting 
that (11+ ± n_)* = a e (%/2<5ai*( 2 ) + *W ± 5oi*(°)) and assuming the normalization J r (*(«))**( m ) = S n , m> the 
gradient terms dependent on 5a± in eq. (|49p are rewritten in the form 



•^GLlgrad = «e£o7 % J dx 



2c: 



m SA x (x) 



V2 



HeSai(x) 



2e 1 1/2 5AJx) - —lm8 ai (x) 

TB 



~j\d x 5ai\ 2 + r~ 2 \8a l \' 



(50) 



In obtaining the | c^Jai | 2 term, the fact that, since J r ($(°))*$( 1 ) = 0, 5ai decouples with 5A X up to the harmonic 
order was used. The last term can be absorbed into the first term of ea. (|4"9"]) so that the mass term vanishes not at T c 
but on the straight H C 2(T) line. As shown in Ref . fl7| , no quadratic terms in 5a\ occur from the sum of the resulting 
|V>| 2 term and fe|-0| 4 /2. 

Here, let us first examine eq.(|50|) by neglecting the x dependences. Then, when 



Sax = 



7 1 l 2 u + i7 1 / 2 u 



V2r B 

where u = uy + vz is the vortex displacement field, the famous Josephson relation 

5A = u x B, 



(51) 



(52) 



or E = —(du/dt) x B implying a nonvanishing vortex flow resistance, follows as a condition for minimizing the 
fluctuation free energy. In fact, this Josephson relation implies the vanishing of the static superfluid rigidity T s j^ to 
a current perpendicular to H, because this relation implies that SA is lost (eaten) by the n = 1 LL fluctuation from 
^GLlgrad so that T s ± oc 8 2 T/5{5A) 2 = 0. 

Next, the x dependences of the fluctuation fields will be incorporated. Further, by substituting eg. ([52")) into the 
magnetic energy term (the last term of eq . (|4"9"|) ~) , the main term B 2 /(4tt) of the vortex tilt modulus C44 is obtained 
from there, while 5C^ defined in sec. I follows from the remaining term in eq. (|50p . i.e., 



al^l\d x 5a!\ 



■a 2 e {d x u) 



B 



1 



B 



(d x u) 2 , 



2r?r ' ' ' 327re 7 (A(T)) 2 V H c2 {T) . 
where A(T) is the London penetration depth. It is easy to see that the resulting expression 

B ( B 



(53) 



1 



167re 7 (A(T)) 2 V H c2 (T) 

is —BM, where M is the magnetization following from the SC condensate. Then, the familiar total tilt modulus 

B(B - 4ttM) BH 



(54) 



An 



-in 



(55) 



is obtained for this Abrikosov vortex lattice. In general, SC44 does not have to be equivalent to —MB: As shown in 
the text, SC44 can become negative in the FFLO vortex lattice, while —M = dT e jdB of the FFLO vortex lattice is 
always positive because the magnetic field tends to destroy superconductivity and hence, increases the condensation 
energy. 



14 



REFERENCES 



[1] A. Bianchi, R.Movshovich, C.Capan, P.G.Pagliuso, and J.L.Sarrao, Phys. Rev. Lett. 91, 187004 (2003). 

[2] C. Capan, A. Bianchi, R. Movshovich, A. D. Christianson, A. Malinowski, M. F. Hundley, A. Lacerda, P. G. Pagliuso, and 

J. L. Sarrao, Phys. Rev. B 70, 134513 (2004). 
[3] C. Martin, C. C. Agosta, S. W. Tozer, H. A. Radovan, E. C. Palm, T. P. Murphy, and J. L. Sarrao, Phys. Rev. B 71, 

020503(R) (2005). 

[4] K. Kakuyanagi, M. Saitoh, K. Kumagai, S. Takashima, M. Nohara, H. Takagi, and Y. Matsuda Phys. Rev. Lett. 94, 
047602 (2005). 

[5] H. Adachi and R. Ikeda, Phys. Rev. B 68, 184510 (2003). 
[6] R. Ikeda and H. Adachi, Phys. Rev. Lett. 95, 269703 (2005). 

[7] P. Fulde and R. A. Ferrell, Phys. Rev. 135, A550 (1964); A. I. Larkin and Y. N. Ovchinnikov, Sov. Phys. JETP 20, 762 
(1965). 

[8] R. Ikeda, |cond-mat/06 10863 (for the Proceedings of M 2 S, Dresden. To appear in Physica C). 



[9] K. Izawa, HYamaguchi, Y.Matsuda, H.Shishido, R.Settai, and Y.Onuki, Phys Rev. Lett. 87, 057002 (2001). 
[10] T. Watanabe, Y. Kasahara, K. Izawa, T. Sakakibara, Y. Matsuda, C. J. van der Beek, T. Hanaguri, H. Shishido, R. Settai, 

and Y. Onuki, Phys. Rev. B 70, 020506(R) (2004). 
[11] This assumption is reasonable in ultrasound measurements with a long enough wavelength of the sound. 
[12] Except in the analysis on the vortex tilt modulus in Appendix B, the magnitude H of the external field H can be identified 

in this paper with the uniform flux density B. 
[13] R.Ikeda and K.Isotani, J. Phys. Soc. Jpn. 67, 983 (1998) ; R. Ikeda, J. Phys. Soc. Jpn. 71, 587 (2002). 
[14] When the vortex lattice with a modulation induced by the paramagnetic depairing has, in contrast to the case of CeCoIns , 

nodal planes parallel to H, the structural transitions between Josephson vortex lattices play important roles even if 

this inequality is satisfied. This case is relevant to the corresponding issue in A-(BET)2FeCU [S. Uji, T. Terashima, M. 

Nishimura, Y. Takahide, T. Konoike, K. Enomoto, H. Cui, H. Kobayashi, A. Kobayashi, H. Tanaka, M. Tokumoto, E. S. 

Choi, T. Tokumoto, D. Graf, and J. S. Brooks, Phys. Rev. Lett. 97, 157001 (2006)] and will be discussed in a separate 

paper. 

[15] PA. Lee and M.G. Payne, Phys. Rev. B 5, 923 (1972). 
[16] R.Ikeda and H. Adachi, Phys. Rev. B 69, 212506 (2004). 

[17] R. Ikeda, T. Ohmi, and T. Tsuneto, J. Phys. Soc. Jpn. 61, 254 (1992). See also R. Ikeda, J. Phys. Soc. Jpn. 64, 3825 
(1995). 

[18] R. Ikeda, |concLmat/0706032"1 l 

[19] M. Houzet and V. P. Mineev, |cond-mat/0703104| 

[20] X. Gratens, L. Mendonca Ferreira, Y. Kopelevich, N. F. Oliveira Jr., P. G. Pagliuso, R. Movshovich, R. R. Urbano, J. L. 
Sarrao, J. D. Thompson, cond-mat/0608722 

[21] We stress that a kink anomaly 2] at -ffFPLO in the thermal conductivity for a heat current || H should appear even when 
the nodal planes are perpendicular to the vortex lines because, in GL region, the thermal conductivity usually increases 
as, just like in entering the FFLO state by increasing the field [lq]. the gap amplitude is reduced. See, for instance, the 
work by C. Caroli and M. Cyrot [ Phys. Kondensierten Materie 4, 285 (1965) ]. Thus, the observation in Ref.2 does not 
contradict the picture with the FFLO state modulating along the field. 

[22] A. I. Larkin and and Y. N. Ovchinnikov, Sov. Phys. JETP 34, 409 (1979). 

[23] H. Won, K. Maki, S. Haas, N. Oeschler, F. Weickert, and P. Gegenwart, Phys. Rev. B 69, 180504(R) (2004). 

[24] See, for instance, U. Klein, D. Rainer, and H. Shimahara, J. Low Temp. Phys. 118, 91 (2000); M. Houzet, A. Buzdin, L. 

Bulaevskii, and M. Maley, Phys. Rev. Lett. 88, 227001 (2002); U. Klein, Phys. Rev. B 69, 134518 (2004); K. Yang and A. 

H. MacDonald, Phys. Rev. B 70, 094512 (2004). 



